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ABSTRACT 


High flux evaporation from a steady meniscus formed in a 2 pm channel is 
modeled using the augmented Young-Laplace equation. The heat flux is found to 
be a function of the long range van der Waals dispersion force which represents 
interfacial conditions between heptane and various substrates. Heat fluxes of (1.3- 
1.6) x 1()6 W/m 2 based on the width of the channel are obtained for heptane 
completely wetting the substrate at 100°C. Small channels are used to obtain these 
large fluxes. Even though the real contact angle is 0°, the apparent contact angle is 
found to vary between 24.8° and 25.6°. The apparent contact angle, which 
represents viscous losses near the contact line, has a large effect on the heat flow 
rate because of its effect on capillary suction and the area of the meniscus. 

The interfacial heat flux is modeled using kinetic theory for the evaporation 
rate. The superheated state depends on the temperature and the pressure of the 
liquid phase. The liquid pressure differs from the pressure of the vapor phase due 
to capillarity and long range van der Waals dispersion forces which are relevant in 
the ultra thin film formed at the leading edge of the meniscus. Important pressure 
gradients in the thin film cause a substantial apparent contact angle for a 
completely wetting system. The temperature of the liquid is related to the 
evaporation rate and to the substrate temperature through the steady heat 
conduction equation. Conduction in the liquid phase is calculated using finite 
element analysis except in the vicinity of the thin film. A lubrication theory 
solution for the thin film is combined with the finite element analysis by the 
method of matched asymptotic expansions. 

KEYWORDS: Micro heat transfer device, thin liquid films, evaporation, 

augmented Young-Laplace equation, contact angle, finite element analysis. 
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INTRODUCTION 


An understanding of the heat transfer characteristics of an evaporating 
extended meniscus is of central importance to the design of many small heat pipes 
and thermosyphons based on interfacial phenomena. Herein, a portion of the 
inverted meniscus system presented in Fig. 1 is analyzed. Using a finite element 
procedure, kinetic theory and the augmented Young-Laplace equation of 
capillarity, the magnitude of the heat flow rate of a steady evaporating "inverted 
meniscus" is numerically evaluated. The long term objective is to determine the 
maximum heat sink capability of the inverted meniscus heat pipe. Since 
resistances associated with fluid flow and evaporation change the shape of the 
evaporating meniscus, the immediate objective is to evaluate the effect of viscous 
stresses on the meniscus profile. Future research will incorporate the results of 
the two-dimensional analysis presented below in a three dimensional study. 

For perspective, the evolution of the present problem can be briefly traced 
as follows. In 1972, Raiff and Wayner [1, 2] experimentally and theoretically 
studied a porous liquid flow control design that was a precursor to the inverted 
meniscus heat pipe studied by Saaski [3]. These studies demonstrated that the 
evaporating inverted meniscus was controlled by the heat transfer characteristics 
of the contact line region. In a complementary theoretical study in 1972, Potash 
and Wayner [4] described an evaporating extended meniscus which consisted of an 
intrinsic meniscus and a thin film extending beyond the apparent contact line of 
the intrinsic meniscus. The two regions are designated in the inverted meniscus 
design presented in Fig. 2. The portion of the extended meniscus that was 
described using the classical equation of capillarity was called the intrinsic 
meniscus. The complementary adjacent thin film was described using the 
disjoining pressure concept of Derjaguin [5, 6]. In 1980, Moosman and Homsy [7] 
presented a procedure which included the effect of conduction in the liquid and 
demonstrated that the relative size of this resistance was substantial. A review ot 
the extensive literature and recent data that leads to the current analysis was 
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recently presented by Wayner and coworkers [8-10]. Recently, Wu and Peterson 
studied a wickless micro heat pipe which had outer dimensions of 1 mm x 1 mm x 
57 mm. Since the fluid had a finite contact angle, they used the classical equation 
of capillarity to describe the internal fluid dynamics [11]. Xu and Carey used a 
completely wetting fluid to study film evaporation from a micro-grooved surface 
with grooves 64 pm wide and 190 pm deep [12]. Although a small evaporating 
meniscus could be a large heat sink, this result has not been conclusively 
demonstrated because of the complex physicochemical phenomena present in the 
microscopic contact line region. 

Existing models apply lubrication theory to the thin film to relate the flow to 
the pressure gradient, apply the augmented Young-Laplace equation to relate the 
pressure in the film to its shape through disjoining pressure and capillarity, and 
link the gradient of flow to the evaporation flux from the film. The flux is linked 
to the local thermodynamic conditions of the liquid phase relative to the vapor 
phase. Use of the augmented Young-Laplace equation is strongly supported by the 
fact that it agrees with the disjoining pressure isotherm for the extreme of a flat 
film and with the Young-Laplace equation for a thick film. Furthermore, it 
provides a convenient mathematical resolution of the breakdown of lubrication 
models based solely on capillary pressure gradients as film thickness approaches 
zero. In the experimental studies of Wayner and coworkers [8, 9] with small 
superheating, the analysis provides a theoretical solution which gives good to 
excellent agreement with the thin flat film data and with the capillary dominated 
thicker (relatively isobaric) portion of the film. There is only one adjustable 
parameter and its behavior is reasonable although its absolute value is not known 
independently. Agreement is more approximate in the transitional portion of the 
films. While the form of the kinetic theory based mass transfer model of 
evaporation is supported, the magnitude of the evaporation coefficient is unknown. 

The extended meniscus in the channel also has a thicker portion in which 
lubrication theory is not valid. In this work, the two regions are analyzed 
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separately and the numerical results are joined together in the framework of the 
method of matched asymptotic expansions to obtain the steady state heat flow rate 
for the complete extended meniscus. Previously, the meniscus in a heat pipe 
evaporator was studied by Stephan and Busse [13]. In their method the meniscus 
shape and heat transfer problem in the micro region and the two dimensional heat 
conduction problem in the macro region were solved. Conduction in the solid 
substrate was also included. Since the solution of each problem was needed as 
input for the other one, an iterative procedure was used. Recent work by 
Swanson and Peterson [14] shed considerable light into the critical mechanisms 
required for proper heat pipe operation. Their thermodynamic model provided 
insights into the behavior of the extended meniscus in capillary structures. It 
included the effects of axial temperature difference, changes in local interface 
curvature, Marangoni effects and disjoining pressure. 

MATHEMATICAL FORMULATION 

Consider a steady liquid meniscus in a two dimensional channel as shown in 
Figs. 1-3. The liquid is superheated relative to the vapor phase by the solid 
surfaces (the slot walls) and so evaporates. A fresh supply is drawn into the 
channel to sustain steady evaporation. Following Wayner and coworkers [8-10] 
the evaporation flux is modeled as 

m = a(T/ v - T v ) + b(P< - P v ) O) 

where a and b depend only on the physical properties of the liquid. The vapor 
phase is assumed to be uniform, having temperature T v and pressure P v . There is 
only one component, both phases are pure. The interfacial temperature of the 
liquid T/v and its pressure P / are in general functions of position. 


The coefficients are 



M y 1/2 | 

i P V M Ah m 'j 

a = C 
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where C is the accommodation coefficient with an ideal value of 2.0, M is the 
molecular weight, R is the universal gas constant, Pv is the bulk vapor pressure of 
the evaporating material at temperature T/ v , V/ is the molar volume of the bulk 
material, Ah m is the enthalpy of vaporization of the material per unit mass at 
temperature T/ V - In Eq.(l), evaporation is promoted by a superheat and hindered 
by a low liquid film pressure. If the temperature differences are small compared 
to the absolute temperature, the known substrate temperature, T s , may be 
substituted for the unknown T/ v in the coefficients. The quantity of interest is the 

heat sink per unit pore depth, Q, 

Q = Ah m } m d l 
C 

where C represents an integration over the length of the liquid-vapor interface. 

The liquid pressure at the interface, P/, is related to the flow of liquid through the 
Navier Stokes and continuity equations. It is also related to the unknown interface 
shape and the vapor phase pressure, P v , through capillarity and long range 
intermolecular forces by the augmented Young-Laplace equation, Eq. (5). 

P/ . p v = - n - a K I 5 ) 


with 


n = 


A 

83 


( 6 ) 


where 6rcA is the Hamaker constant (negative for a completely wetting liquid), A 
is the dispersion constant, 5(x) is the film thickness, K is the curvature and a is the 
surface tension. The first term on the rhs of Eq. (5) is called the disjoining 
pressure, n, and it represents the change in the body force on the liquid due to the 
long range van der Waals forces between the liquid and solid over a narrow range 
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of thicknesses. Herein, a "pressure jump" model is used to account for both the 
disjoining pressure and capillarity. The disjoining pressure is really a model of the 
film thickness dependent intermolecular force field. The vapor is assumed to be 
pure. These forces are relevant in the very thin portions of the liquid phase. The 
temperature T/ v is related to the constant wall temperature T s through the steady 
diffusion equation, convective effects taken to be negligible. The liquid is assumed 
to completely wet the pure wall; forming a microscopic contact angle of 0°. 
Therefore, a portion of the liquid phase may be analyzed with lubrication theory 
due to the small slope. In this portion the natural length scale is the thickness of 
the thin film whereas in the rest of the liquid phase the natural length scale is the 
width of the channel, W. The radical divergence of these length scales suggests 
analysis with matched asymptotic expansions. 


Lubrication Region, Q| u b 


The literature provides relationships which are the building blocks for the 
following lubrication model of the evaporating thin film portion [see e.g., 8-10]. 
The liquid pressure at the interface, P / , is related to the vapor phase pressure, P v , 
by Eqs. (5) and (6). The curvature is simply the second derivative of the film 
thickness with distance along the substrate. Lubrication theory for a slightly 
tapered film relates the mass flow rate in the thin film to the pressure gradient in 
the direction of flow: 


r(x) 


-83 

3v 



(7) 


where T(x) is the mass flow rate and v is the kinematic viscosity. These two 
equations demonstrate that the mass flow rate, T(x), is a function of the thin film 
shape, 5(x). The evaporative flux, m, is modeled according to Eq. (1). 

Following Moosman and Homsy [7], the one dimensional conduction heat 
transfer solution for the film is used to eliminate T/v in favor of Ts: 
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m = — [a(T s - T v ) + b( -n - aK )] (8) 

aAhm „ 

1 + nr -6 

where k is the thermal conductivity of the liquid. This equation clearly 
demonstrates the substantial effect of film thickness on the evaporation rate. 

Therefore, both fluid flow and evaporation depend on the same 
intermolecular force field which is a function of the film profile. 


Using a steady state material balance, the evaporative flux is 

ar 


(9) 


The coupled differential equations are found to be Eq. (5) and the following 
combination of Eqs. (7), (8), and (9): 


3v d x ^ dx J 


1 + 


_J 

aAhm 

k 


5 


[aAT + b( -n - aK )] 


( 10 ) 


where AT = T s - T v . 


The above equations are satisfied by the interline film thickness: 

M / 3 


So — 


- A b 
a (T s - T v ) 


(ID 


which provides the farfield condition 
5 — > 8 0 , x -> -°° 

These equations have been previously combined and solved numerically for 
the film thickness profile and the evaporation rate, T(x) by Wayner and coworkers 
[8,9]. In this work, the Adams method solver was used to find a solution satisfying 
a second farfield condition 

~ tan 9 e x — > (12) 

d x 
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where the effective or macroscopic apparent contact angle, 0 e , is thus defined. 

This quantity depends only on the physical properties of the liquid, the dispersion 
constant. A, and the temperature difference (T s - T v ). We note that the apparent 
contact angle, 0 e , for this completely wetting system represents the viscous losses 
in the extremely thin portion of the meniscus. The contribution from the thin film 
lubrication region, Q| u b> to the total heat sink, Q, was calculated from T(x). The 
solution for the profile , 5(x), which is schematically presented in Fig. 4, is the 
leading order approximation of the interface shape in the lubrication region. 

Macroscopic Region, Qintr 

The characteristic length of the macroscopic region is the channel width W. 
The characteristic fluid velocity may be estimated from Eq.(l) as 

V= a (T s - T v ) / p (13) 

which is an upper bound. Using a viscous scaling appropriate for a small channel, 
the pressure change over a length W due to this flow is scaled as 

AP = a (T s - T v ) v / W (14) 

The capillary pressure jump at the interface is scaled as 

AP C = a / (W/2) (15) 

For the case of interest the viscous pressure scale, AP, is much smaller than the 
capillary pressure scale, AP C , by three orders of magnitude. Therefore the changes 
in the liquid pressure P / from place to place on the liquid vapor interface due to 
viscous flow in the macroscopic region are utterly insignificant in comparison to 
the capillary pressure jump. Therefore, the leading order approximation of the 
curvature of the liquid vapor interface must be a constant. Thus, the interface is 
described by a section of a circle. This interface has an asymptotic form valid near 
a channel wall which is illustrated in Fig. 5. 
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Consider the circular arc: 


, \ i 

( x - x 0 ) + 



( 16 ) 


centered around a point located midway between the channel walls as shown in 
Fig. 5. The walls are described by the planes y = 0 and y = W. This arc is analyzed 
near y = 0 with a Taylor expansion to yield 


y 




+ O (x 2 /W) 


x — > 0 


(17) 


Since the lubrication equation is translationally invariant with respect to x, the 
farfield form of the lubrication solution, Equation (12), can be written as 

8iub —* x tan 0 e x— (18) 


Therefore, if R is chosen to be 

W 1 

R - Q 

2 cos 9 e 


(19) 


the outer limit of the lubrication solution matches the inner limit of the 
macroscopic solution. Therefore, the shape of the meniscus is defined by the pore 
geometry and the effective contact angle 0 e- 


The contribution of the macroscopic solution in the intrinsic region, Qintr- t0 
the total heat sink, Q, is found by solving 


V 2 T = 0 
T = T S 


( 20 ) 


on the channel walls 

on the truncation segments and on a bounding 
surface through the liquid phase, well removed from the meniscus. 


3T 

-k — = 0 
dn 


3T a on the liquid vapor interface where the 

3 ^ = a (T ‘ Tv) - b ^ outward unit normal vector is used. 
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There is a singularity where the circular arc intersects the channel wall which can 
be avoided by using the lubrication solution. The numerical procedures to 
accomplish this are discussed below. 

The steady state heat equation was solved by finite element analysis. The 
singularity was handled by truncating the macroscopic region at a finite film 
thickness as shown in Fig. 6b. Since heat transfer in the lubrication region is one 
dimensional and perpendicular to the channel walls, the truncation surfaces were 
taken to be insulated. Even though the slot is semi-infinite in the "x" direction, 
heat transfer was anticipated only within two or three channel widths of the 
interface. Therefore, the domain was terminated with an insulated surface which 
appears on the right hand side of Fig 6a. Geometric modeling of the boundaries 
and mesh generation were performed with software developed at the Scientific 
Computation Research Center at Rensselaer. One of the meshes used is shown in 
Figs. 6a and b with magnification of the smaller scale portions of the liquid phase. 
Finite element analysis was performed with ABAQUS and a user defined element. 
The effect of pressure on m was estimated to be negligible in the macroscopic 
region. The macroscopic contribution to the heat sink Q was obtained by 
integration using the values of T at the mesh nodes on the boundary. 

The dimensionless interfacial temperature is displayed as a function of the 
dimensionless distance along the channel wall in Fig. 7. A dimensionless 
temperature of 1.0 corresponds to the substrate temperature T s at a film thickness 
of 5 0 , whereas a dimensionless value of 0 corresponds to the vapor temperature T v . 
In the macroscopic region, the interfacial mass flux can be obtained using these 
results with Eq. (1) while neglecting the effect of pressure. The lubrication region 
ends where the film thickness is roughly 20-30 times 8 0 . The scaling of distance is 
shown in Fig. 6a where a dimensional distance of 2 pm is equal to a scaled distance 
of 4. The macroscopic region begins at -5.968. Note the mesh and model in Figs. 

6a and b corresponds to the results displayed in Fig. 7. 
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Two meshes were tested in a sample calculation: a fine mesh and a coarse 
mesh. The variation in the line heat sink was only 0.1%. The element edges in the 
fine mesh are approximately 60% of the length of the element edges in the coarse 
mesh. Two cutoffs were also studied in a sample computation. The cutoffs were at 
+6 and +3 units (defined by the coordinate system displayed in Fig. 6a). The 
difference in the macroscopic line sink was 2 (10'3)%. 

DISCUSSION 

The results for various values of the modified Hamaker constant are given in 
Table 1. The range of Hamaker values represent a clean silicon surface at the 
higher value and a contaminated surface at the lower value. The intermediate 
value represents Si02- The large vapor pressure value is selected to enhance the 
heat flux. The predicted high average heat fluxes (based on the channel width) are 
apparently the result of the high vapor pressure, the large superheat, and the 
small channel width. The "thicker" portions of the interface are relatively inactive 
as shown by Figs. 7 and 6a. This is due to the assumption or rather design of a 
constant temperature at the channel walls and to the ease of conduction through 
the thinner portion of the liquid region relative to the thicker portions. Therefore, 
a larger channel with the same temperature difference, T s - T v , would probably 
have a reduced average evaporative heat flux. 

The algorithm is attractive because it is not necessary to re-solve the 
lubrication portion if the channel width is changed, provided the basic assumptions 
are met (i.e., the relative width of the channel is sufficiently large.) Furthermore, an 
iterative computational scheme joining the finite element analysis and an ordinary 
differential equation solver is avoided. Both the finite element code and the 
differential equation solver used are standard. Finally the macroscopic apparent 
contact angle, 0e, is identified as a fairly fundamental quantity of great importance 
affecting both the area of the meniscus and of course the capillary suction of the 
device. Although the variation of 0 e in these case studies is smaller than the 



variation of the dispersion constant, that constant is susceptible to wide variations 
as shown by the fact that not all fluids completely wet a given substrate. 

PRACTICAL SIGNIFICANCE AND FUTURE RESEARCH 

The present work has theoretically studied the high heat transfer potential 
of a micro heat pipe and has showed that a heat flux of 1 .3-1.6 x 10 6 W/m 2 is 
possible from a steady meniscus in a small channel. It demonstrates the 
importance that interline phenomena play in the emerging area of micro heat 
transfer. The model and the information obtained in this theoretical study have 
enhanced the understanding of the complex physicochemical phenomena 
associated with high heat flux evaporation processes and can be of importance in 
designing practical micro heat transfer devices. For example, since the change in 
the apparent contact angle due to the viscous stress in the contact line region 
reduces the suction potential at the base of the meniscus, the evaporating 
meniscus becomes unstable at a lower heat flux than would be predicted without 
the resistance. 

Even though the calculations suggest that very high heat flux may be 
sustained by an evaporating meniscus in a small region, scaling up to a three 
dimensional device with a myriad of channels is not a trivial exercise. This is due 
to the fact that the pressure drop associated with the vapor removal would be 
substantial and a function of design. Furthermore, we note that it is 
experimentally very difficult to optically study the evaporating meniscus for such 
a high degree of superheat and the whole issue of stability at high heat flux at 
extremely small systems is yet to be addressed and should be considered in future 
studies. There are many unknowns in this small region. 

CONCLUSIONS 

1) A stable evaporating meniscus with a very high heat flux based on the width 

of the channel is theoretically possible. 



2) The simulation suggests that very high heat fluxes based on the width of the 
channel are associated with micron sized channels. 

3) The extremely large theoretical heat fluxes in the lubrication region have not 
been studied experimentally. 

4) A large variation in the Hamaker constant was not associated with an equally 
large variation in the total heat sink or the apparent contact angle. 

5) The large apparent contact angle for a completely spreading system is 
(probably) due to large viscous stresses in the contact line region where fluid 
flows in an extremely thin film. 

6) The apparent contact angle has a large effect on the area of the meniscus and 
the capillary suction of the device. 
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NOMENCLATURE 

A Hamaker constant, J 

A dispersion constant, J 

a defined by Eq. (2), kg/m-Ks 

b defined by Eq. (3), s/m 

C evaporation coefficient, dimensionless 



Ah m molar heat of vaporization, J/mol 
k thermal conductivity, W/mK 

K curvature, nr 1 

M molecular weight, g/mol 

m interfacial mass flux, kg/m 2 s 

n normal direction, m 

P pressure, N/m 2 

Q heat flow rate, W/m 

R gas constant, J/Kmol 

R radius of curvature, m 

T temperature, K 

V characteristic fluid velocity, m/s 

V molar volume, m 3 /mol 

W channel width, m 

x coordinate distance, m 

y coordinate distance, m 

8 film thickness, m 

A difference 

T mass flow rate, kg/m 

v kinematic viscosity, m 2 /s 

a surface tension, N/m 

n disjoining pressure, N/m 2 

p density, kg/m 3 

0 contact angle, deg 

Subscripts 

c capillary 

E apparent 

intr intrinsic region 

/ liquid 

lub lubrication region 
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o 


contact line 

/v liquid-vapor interface 

s solid 

Total total of lub+intr 
v vapor 
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TABLE 1 


SUMMARY OF HEAT PIPE CASE STUDIES 


Some of the Inputs 

-A (J) 

T (K) 

AT (K) 
v (m 2 /sec) 
ct (N/m) 


Case: 1 Case: 

10-21 10-22 

373 373 

5.19 5.19 

3.233 x 10- 7 3.233 

1.247 x 10-2 1.247 


Case: 3 

10-23 

373 

5.19 

10- 7 3.233 x lO* 7 

10-2 1.247 x lO* 2 


Calculated Values 


a (kg/m 2 »K*sec) 


b (sec/m) 

§0 (A) 


13.19 

2.536 x lO- 5 
7.18 


mjd = aAT (kg/m 2 »sec) 68.456 


13.19 

2.536 x 10-5 
3.33 

68.456 


13.19 

2.536 x 10-5 
1.55 

68.456 


Outputs 


0E 

Qiub. (W/m) 

Qintr. (W/m) 

QTotal (W/m) 

HEAT FLUX (W/m 2 ) 
based on 
W = 2 x 10- 6 m 


23.75 

1.238 

1.883 

3.121 

1.5605 x I0 6 


24.80 
0.7396 
1.986 
2.726 _ 

1.363 x 10 6 


25.43 

0.460 

2.087 

2.547 

1.2735 x 106 
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Fig. 1. 
Fig. 2. 


Fig. 3. 
Fig. 4. 

Fig. 5. 
Fig. 6a. 
Fig. 6b. 
Fig .7 


Inverted meniscus heat pipe concept. 

Schematic of Inverted Extended Evaporating Meniscus. Liquid enters 
through the porous flow control device that prevents the vapor from 
entering the liquid space above the flow control device. The passive 
shape of the liquid meniscus, which is a function of the transport 
processes and interfacial phenomena, responds to the pressure 
distribution in the liquid and vapor spaces which, in turn, responds to 
the heat flux distribution. The interfacial physics of the intermolecular 
force field in the extended meniscus which form the basis for the 
operation of this device are discussed in [4,10]. The vapor exits in the 
direction perpendicular to the plane of the paper. 

Steady state evaporating meniscus in a channel. 

Schematic diagram of solution on the scale of the adsorbed thin film 
thickness, 6 0 . The apparent contact angle, ©e, is a function of the 
properties of the liquid, the dispersion constant, and the superheat, 

T s - T v . 

Solution on the scale of the channel width, W. Circular arc centered at 
(x 0 , W/2) . The radius R is chosen to provide the proper slope at the 
walls using Eq. (17). 

Finite element mesh (fine mesh) for solution of steady heat equation on 
macroscopic scale, each grid block (dotted lines) is 0.5 micron on edge. 
Detail of leading edge is shown in Figure 6b. 

Finite element mesh (fine mesh) at leading edge of the macroscopic 
region. The vertical edge is the truncation surface, 7.495 nm high, the 
rest of the interface is modeled with lubrication theory. 

Dimensionless temperature, (T]v - T V )/(T S - T v ), versus dimensionless 
distance, defined in Figure 6a. 



Fig 1 



V : VAPOR FLOW CHANNELS . L : LIQUID , S : SOLID FIN 
H CRITICAL SHAPE CONTROLS THE HEAT FLUX 

Land stability of the process 










Fig 2 


LIQUID FLOW 



POROUS 


LIQUID 


SOLID 



VAPOR 


HSBBOsrA 




THICKER INTRINSIC MENISCUS REGION 
CONTROLLED BY CAPILLARITY 


THIN FILM REGION CONTROLLED 
BY DISJOINING PRESSURE AND 
CAPILLARITY 


















